Última actualización: 20 mayo 2020
En esta entrega se hace una verificación complementaria en función de:
La siguiente figura muestra un resumen de la producción anual de sedimentos ordenado según el uso de la tierra. Los usos EUCA y MONT tienen una producción de sedimentos bastante inferior AGRP, AGRC, LECH y GRAS. Adicionalmente AGRP, LECH, y GRAS tienen un comportamiento bastante similar. AGRC está ligeramente por encima y su variabilidad practicamente duplica la variabilidad del grupo AGRP, LECH y GRAS.
SYLD <- readRDS("~/R_projects/SWAT_calibration_ETmodis/resultados/optim20200401_22/SYLD_promedio_10anos.RDS")
hru_info <- readRDS("~/R_projects/SWAT_calibration_ETmodis/shapefiles_RDS/hru_info.RDS")
SYLD_GRAS = unlist(SYLD[(hru_info$LANDUSE=="GRAS")])
SYLD_EUCA = unlist(SYLD[(hru_info$LANDUSE=="EUCA")])
SYLD_MONT = unlist(SYLD[(hru_info$LANDUSE=="MONT")])
SYLD_AGRC = unlist(SYLD[(hru_info$LANDUSE=="AGRC")])
SYLD_AGRP = unlist(SYLD[(hru_info$LANDUSE=="AGRP")])
SYLD_LECH = unlist(SYLD[(hru_info$LANDUSE=="LECH")])
library(plotly)
fig <- plot_ly(y = ~SYLD_AGRC, type = "box", name="AGRC")
fig <- fig %>% add_trace(y = ~SYLD_AGRP, name="AGRP")
fig <- fig %>% add_trace(y = ~SYLD_LECH, name="LECH")
fig <- fig %>% add_trace(y = ~SYLD_GRAS, name="GRAS")
fig <- fig %>% add_trace(y = ~SYLD_EUCA, name="EUCA")
fig <- fig %>% add_trace(y = ~SYLD_MONT, name="MONT")
fig <- fig %>% layout(yaxis = list(title="Producción anual de sedimentos [ton/ha]"),
showlegend = FALSE,
autosize=FALSE,
width=500,
height=350)
figA continuación se muestra la producción de sedimentos en funcion de las 6 clases de pendientes modeladas en SWAT:
## [1] "S1(0-3.0)" "S2(3.0-6.0)" "S3(6.0-9.0)" "S4(9.0-12.0)"
## [5] "S5(12.0-15.0)" "S6(15.0-9999)"
slope_plot = function(Luse){
c1 = hru_info$LANDUSE==Luse & hru_info$SLOPE_BAND==Sband[1]
S1 = unlist(SYLD[c1])
c2 = hru_info$LANDUSE==Luse & hru_info$SLOPE_BAND==Sband[2]
S2 = unlist(SYLD[c2])
c3 = hru_info$LANDUSE==Luse & hru_info$SLOPE_BAND==Sband[3]
S3 = unlist(SYLD[c3])
c4 = hru_info$LANDUSE==Luse & hru_info$SLOPE_BAND==Sband[4]
S4 = unlist(SYLD[c4])
c5 = hru_info$LANDUSE==Luse & hru_info$SLOPE_BAND==Sband[5]
S5 = unlist(SYLD[c5])
c6 = hru_info$LANDUSE==Luse & hru_info$SLOPE_BAND==Sband[6]
S6 = unlist(SYLD[c6])
if(is.null(S1)){S1=0}
if(is.null(S2)){S2=0}
if(is.null(S3)){S3=0}
if(is.null(S4)){S4=0}
if(is.null(S5)){S5=0}
if(is.null(S6)){S6=0}
fig <- plot_ly(y = ~S1, type = "box", name=paste0("s1 (",table(c1)["TRUE"],")" ))
fig <- fig %>% add_trace(y = ~S2, name=paste0("s2 (",table(c2)["TRUE"],")" ))
fig <- fig %>% add_trace(y = ~S3, name=paste0("s3 (",table(c3)["TRUE"],")" ))
fig <- fig %>% add_trace(y = ~S4, name=paste0("s4 (",table(c4)["TRUE"],")" ))
fig <- fig %>% add_trace(y = ~S5, name=paste0("s5 (",table(c5)["TRUE"],")" ))
fig <- fig %>% add_trace(y = ~S6, name=paste0("s6 (",table(c6)["TRUE"],")" ))
fig <- fig %>% layout(title=paste(Luse),
yaxis = list(title="Producción anual de sedimentos [ton/ha]",
range = c(0, 70)),
showlegend = FALSE,
autosize=FALSE,
width=350,
height=300)
fig
}
leafsync::sync(slope_plot("AGRC"),slope_plot("AGRP"),
slope_plot("LECH"),slope_plot("GRAS"))La verificación de cantidad de agua se hace a escala mensual en el período 2006-2015 (2005 de calentamiento). Esta verificacion muetra a grandes trazos que en los primeros años de simulacion (2006-2010) predomina la subestimación. Luego a partir del año 2010 se sobreestima el caudal.
A continuación leemos los resultados del modelo ejecutado en ese período así como los valores observados. Nótese que la simulación de caudal tiene de 100 ejecuciones del modelo, los estadísticos se dan en funcion de la mediana de las 100 simulaciones.
flow_m23 <- readRDS("~/R_projects/SWAT_calibration_ETmodis/resultados/optim20200401_22/flow_m23.RDS")
flow_m30 <- readRDS("~/R_projects/SWAT_calibration_ETmodis/resultados/optim20200401_22/flow_m30.RDS")
flow_m36 <- readRDS("~/R_projects/SWAT_calibration_ETmodis/resultados/optim20200401_22/flow_m36.RDS")
flow_m38 <- readRDS("~/R_projects/SWAT_calibration_ETmodis/resultados/optim20200401_22/flow_m38.RDS")
flow_m41 <- readRDS("~/R_projects/SWAT_calibration_ETmodis/resultados/optim20200401_22/flow_m41.RDS")
flow_m23 = lapply(flow_m23[,-1],function(x){convertFlow(zoo(x, flow_m23$date),
from="cumecs", to="mm", area.km2 = 687,
timestep.default = "months")})
flow_m30 = lapply(flow_m30[,-1],function(x){convertFlow(zoo(x, flow_m30$date),
from="cumecs", to="mm", area.km2 = 1092,
timestep.default = "months")})
flow_m36 = lapply(flow_m36[,-1],function(x){convertFlow(zoo(x, flow_m36$date),
from="cumecs", to="mm", area.km2 = 2744,
timestep.default = "months")})
flow_m38 = lapply(flow_m38[,-1],function(x){convertFlow(zoo(x, flow_m38$date),
from="cumecs", to="mm", area.km2 = 3159,
timestep.default = "months")})
flow_m41 = lapply(flow_m41[,-1],function(x){convertFlow(zoo(x, flow_m41$date),
from="cumecs", to="mm", area.km2 = 4896,
timestep.default = "months")})
Qsim = cbind(Qsim_23 = apply(as.data.frame(flow_m23),1,median),
Qsim_30 = apply(as.data.frame(flow_m30),1,median),
Qsim_36 = apply(as.data.frame(flow_m36),1,median),
Qsim_38 = apply(as.data.frame(flow_m38),1,median),
Qsim_41 = apply(as.data.frame(flow_m41),1,median))
Qobs <- readRDS("~/R_projects/SWAT_calibration_ETmodis/data/runoff_mensual_2006_2015.RDS")
qo = cbind(Qobs, Qsim)La siguiente tabla muestra los estadísticos KGE, NSE y Pbias a escala mensual en el período 2006-2015.
library(hydroGOF)
kge_stat = round(c(KGE(qo$Qsim_23,qo$Qobs_23),
KGE(qo$Qsim_30,qo$Qobs_30),
KGE(qo$Qsim_36,qo$Qobs_36),
KGE(qo$Qsim_38,qo$Qobs_38),
KGE(qo$Qsim_41,qo$Qobs_41)),2)
nse_stat = round(c(NSE(qo$Qsim_23,qo$Qobs_23),
NSE(qo$Qsim_30,qo$Qobs_30),
NSE(qo$Qsim_36,qo$Qobs_36),
NSE(qo$Qsim_38,qo$Qobs_38),
NSE(qo$Qsim_41,qo$Qobs_41)),2)
pbias_stat = round(c(pbias(qo$Qsim_23,qo$Qobs_23),
pbias(qo$Qsim_30,qo$Qobs_30),
pbias(qo$Qsim_36,qo$Qobs_36),
pbias(qo$Qsim_38,qo$Qobs_38),
pbias(qo$Qsim_41,qo$Qobs_41)),1)
stat = rbind(KGE=kge_stat,NSE=nse_stat, PBIAS=pbias_stat)
kable(stat, col.names = c("WL_23", "WL_30", "WL_36", "WL_38", "WL_41")) %>%
kable_styling(c("striped", "bordered"))| WL_23 | WL_30 | WL_36 | WL_38 | WL_41 | |
|---|---|---|---|---|---|
| KGE | 0.25 | 0.32 | 0.74 | 0.37 | 0.70 |
| NSE | -0.26 | 0.34 | 0.73 | 0.33 | 0.79 |
| PBIAS | 47.20 | -16.10 | 21.60 | -46.50 | 26.90 |
res = cbind(r23 = (qo$Qsim_23 - qo$Qobs_23),
r30 = (qo$Qsim_30 - qo$Qobs_30),
r36 = (qo$Qsim_36 - qo$Qobs_36),
r38 = (qo$Qsim_38 - qo$Qobs_38),
r41 = (qo$Qsim_41 - qo$Qobs_41))
Res_plot = function(x){
cc = substr(index(res),1,3) %in% "ene"
nn = index(res)
nn[!cc] = NA
nn = substr(nn,6,10)
barplot(res[,x], las=1, col=ifelse(res[,x]>=0,4,2),
main = x, ylab="Qsim-Qobs (mm)", names=nn, ylim = c(-100,100))
}
par(mfrow=c(5,1), mar=c(4,4,1,1))
Res_plot("r23")
Res_plot("r30")
Res_plot("r36")
Res_plot("r38")
Res_plot("r41")A work by Rafael Navas